using System;
using System.Linq;

using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Physics;
using Atomic.Structures;
using Atomic.Thermodynamics.StateEquations;

namespace Atomic.Thermodynamics
{
	[Serializable]
	public class DebyeGruneisenPotential : IHelmholtzPotential
	{
		private double b0, v0, m, scaling, gamma, vs, theta0;
		private int n;

		public DebyeGruneisenPotential(double equilibriumBulkModulus, double equilibriumVolume, double mass, int atoms, double scaling, double gamma)
		{
			if (equilibriumBulkModulus <= 0.0 || equilibriumVolume <= 0.0 || mass <= 0.0 || atoms < 1 || scaling <= 0.0 || gamma < 0.0)
			{
				throw new ArgumentOutOfRangeException();
			}

			// Debye-Grüneisen model scaling parameter.
			// 0.617 recommended by Moruzzi, PRB 37 (1988) 790.
			//double scaling = 0.617;

			// 1 for low temperature, 2/3 for high temperature.
			//double gamma0 = 2.0 / 3.0;

			// Gruneisen constant.
			//double gamma = (1.0 + equilibriumBulkModulusVolumeDerivative) / 2.0 - gamma0;

			b0 = equilibriumBulkModulus;
			v0 = equilibriumVolume;
			m = mass;
			n = atoms;
			this.scaling = scaling;
			this.gamma = gamma;

			// The "natural" units of mass in this framework is eV*s^2/Å^2 (cf. 6).
			// Obtain mass by converting to SI first (1 u = 1.66e-27 kg, 1 eV = 1.602e-19 J).
			//double m = mass * 1.66053892173e-27 / (1.60217656535e-19 / 1.0e-20);

			double hbar = 6.5821192815e-16; // eV*s
			double kb = 8.617332478e-5; // eV/K

			// Sound velocity (Å/s).
			vs = scaling * Math.Sqrt(b0 * v0 / m);

			// Debye temperature (K).
			theta0 = hbar / kb * Math.Pow(6.0 * Math.PI * Math.PI * atoms / v0, 1.0 / 3.0) * vs;
		}

		public static double StructureMass(Structure structure)
		{
			return structure.RemoveVacancies().Sites.Sum(s => s.Atom.Mass) * Conversions.AtomicMassUnit;
		}

		public static double StructureMassGeometric(Structure structure)
		{
			structure = structure.RemoveVacancies();

			double atoms = structure.Sites.Count;
			double geometricAverage = Math.Exp(structure.Sites.Sum(s => Math.Log(s.Atom.Mass)) / atoms);
			return atoms * geometricAverage * Conversions.AtomicMassUnit;
		}

		public double FreeEnergy(double temperature, double volume)
		{
			double f, f_v, f_vv, f_vt, s, cv;
			Compute(temperature, volume, out f, out f_v, out f_vv, out f_vt, out s, out cv);

			return f;
		}

		public double Pressure(double temperature, double volume)
		{
			double f, f_v, f_vv, f_vt, s, cv;
			Compute(temperature, volume, out f, out f_v, out f_vv, out f_vt, out s, out cv);

			// http://en.wikipedia.org/w/index.php?title=Relations_between_heat_capacities&oldid=579503942#Derivation
			return -f_v;
		}

		public double PressureTemperatureDerivative(double temperature, double volume)
		{
			double f, f_v, f_vv, f_vt, s, cv;
			Compute(temperature, volume, out f, out f_v, out f_vv, out f_vt, out s, out cv);

			return -f_vt;
		}

		public double BulkModulus(double temperature, double volume)
		{
			double f, f_v, f_vv, f_vt, s, cv;
			Compute(temperature, volume, out f, out f_v, out f_vv, out f_vt, out s, out cv);

			return volume * f_vv;
		}

		/// <summary>
		/// The Debye temperature (K) given volume (Å^3) as corrected by the Debye-Grüneisen model.
		/// </summary>
		public double DebyeTemperature(double volume)
		{
			if (volume <= 0.0)
			{
				throw new ArgumentOutOfRangeException();
			}

			return theta0 * Math.Pow(v0 / volume, gamma);
		}

		private void Compute(double t, double v, out double f, out double f_v, out double f_vv, out double f_vt, out double s, out double cv)
		{
			if (t < 0.0 || v <= 0.0)
			{
				throw new ArgumentOutOfRangeException();
			}

			//double theta0 = scaling * 231.0394 * Math.Pow(v0, 1.0 / 6.0) * Math.Sqrt(b / m);
			//double theta0 = DebyeTemperature;

			// Debye temperature at given volume in the Debye-Grüneisen model.
			double theta = theta0 * Math.Pow(v0 / v, gamma);
			double theta_v = -gamma * theta / v;
			double theta_vv = -gamma * (theta_v / v - theta / (v * v));
			//double theta_vt = 0.0;

			double x = theta / t;
			double x_v = theta_v / t;
			double x_vv = theta_vv / t;
			double x_t = -theta / (t * t);
			double x_vt = -theta_v / (t * t);
		
			// The Debye function D_3 and its derivatives.
			double d, d_v, d_vv, d_vt;
			DebyeFunction(x, x_v, x_vv, x_t, x_vt, out d, out d_v, out d_vv, out d_vt);

			if (t == 0.0)
			{
				// The last factor dominates as the product converges towards 0. Guessed by numerial testing.
				d_v = 0.0;
				d_vv = 0.0;
				d_vt = 0.0;
			}

			// Boltzmann's constant.
			double kb = 8.617332478e-5; // eV/K

			double e = kb * (3.0 * t * d + 1.125 * theta);
			double e_v = kb * (3.0 * t * d_v + 1.125 * theta_v);
			double e_vv = kb * (3.0 * t * d_vv + 1.125 * theta_vv);
			double e_vt = kb * 3.0 * (1.0 * d_v + t * d_vt);

			double z = 1.0 - Math.Exp(-x);
			double z_v = (1.0 - z) * x_v;
			double z_vv = (z - 1.0) * x_v * x_v + (1.0 - z) * x_vv;
			double z_t = (1.0 - z) * x_t;
			double z_vt = (z - 1.0) * x_t * x_v + (1.0 - z) * x_vt;

			s = kb * (4.0 * d - 3.0 * Math.Log(z));
			double s_v, s_vv, s_vt;
			if (t == 0.0)
			{
				// Limiting values guessed by numerial testing.
				s_v = 0.0;
				s_vv = 0.0;
				s_vt = 0.0;
			}
			else
			{
				s_v = kb * (4.0 * d_v - 3.0 / z * z_v);
				s_vv = kb * (4.0 * d_vv + 3.0 / (z * z) * z_v * z_v - 3.0 / z * z_vv);
				s_vt = kb * (4.0 * d_vt + 3.0 / (z * z) * z_t * z_v - 3.0 / z * z_vt);
			}

			f = e - t * s;
			f_v = e_v - t * s_v;
			f_vv = e_vv - t * s_vv;
			f_vt = e_vt - (1.0 * s_v + t * s_vt);

			cv = 3.0 * kb * (4.0 * d - 3.0 * x / (Math.Exp(x) - 1.0));
		
			// Perform scaling.
			f *= n;
			f_v *= n;
			f_vv *= n;
			f_vt *= n;
			s *= n;
			s_v *= n;
			s_vv *= n;
			s_vt *= n;
			cv *= n;
		}

		public static double DebyeFunction(double x)
		{
			if (x < 0.0)
			{
				return double.NaN;
			}

			if (x == 0.0)
			{
				// http://en.wikipedia.org/w/index.php?title=Debye_function&oldid=548162034#Limiting_values
				return 1.0;
			}

			if (double.IsPositiveInfinity(x))
			{
				return 0.0;
			}

			return Gsl.gsl_sf_debye_3(x);
		}

		public static double DebyeFunctionDerivative(double x)
		{
			//return 3.0 * (-DebyeFunction(x) / x + 1.0 / (Math.Exp(x) - 1.0));
			return DebyeFunctionDerivative(x, DebyeFunction(x));
		}

		public static double DebyeFunctionDerivative(double x, double d)
		{
			return 3.0 * (-d / x + 1.0 / (Math.Exp(x) - 1.0));
		}

		private static void DebyeFunction(double x, double x_v, double x_vv, double x_t, double x_vt, out double d, out double d_v, out double d_vv, out double d_vt)
		{
			double z = Math.Exp(x) - 1.0;
			double z_v = (z + 1.0) * x_v;
			//double z_vv = z_v * x_v + (z + 1.0) * x_vv;
			double z_t = (z + 1.0) * x_t;
			//double z_vt = z_t * x_v + (z + 1.0) * x_vt;

			d = DebyeFunction(x);
			d_v = 3.0 * (-d / x + 1.0 / z) * x_v;
			double d_t = 3.0 * (-d / x + 1.0 / z) * x_t;
			d_vv = 3.0 * (-d_v / x + d / (x * x) * x_v - 1.0 / (z * z) * z_v) * x_v + 3.0 * (-d / x + 1.0 / z) * x_vv;
			d_vt = 3.0 * (-d_t / x + d / (x * x) * x_t - 1.0 / (z * z) * z_t) * x_v + 3.0 * (-d / x + 1.0 / z) * x_vt;
		}

		/// <summary>
		/// Equilibrium bulk modulus (eV/Å^3).
		/// </summary>
		public double EquilibriumBulkModulus
		{
			get
			{
				return b0;
			}
		}

		/// <summary>
		/// Equilibrium volume of the unit cell (Å^3).
		/// </summary>
		public double EquilibriumVolume
		{
			get
			{
				return v0;
			}
		}

		/// <summary>
		/// Total mass of atoms (ions) in the unit cell (eV*s^2/Å^2, 1 u = ???).
		/// </summary>
		public double Mass
		{
			get
			{
				return m;
			}
		}

		/// <summary>
		/// Number of atoms (ions) in the unit cell.
		/// </summary>
		public int Atoms
		{
			get
			{
				return n;
			}
		}

		/// <summary>
		/// Scaling parameter of sound velocity.
		/// </summary>
		public double Scaling
		{
			get
			{
				return scaling;
			}
		}

		/// <summary>
		/// Gruneisen constant.
		/// </summary>
		public double Gamma
		{
			get
			{
				return gamma;
			}
		}

		/// <summary>
		/// Computed sound velocity (Å/s).
		/// </summary>
		public double SoundVelocity
		{
			get
			{
				return vs;
			}
		}

		/// <summary>
		/// The computed Debye temperature at equilibrium volume.
		/// </summary>
		public double EquilibriumDebyeTemperature
		{
			get
			{
				return theta0;
			}
		}

		IGibbsPotential IHelmholtzPotential.GibbsPotential
		{
			get
			{
				// Is it useful to implement this at all?
				throw new NotImplementedException();
			}
		}

		double IHelmholtzPotential.EquilibriumVolume(double temperature, double pressure)
		{
			if (pressure != 0.0)
			{
				throw new NotImplementedException();
			}

			// Not just EquilibriumVolume??
			//return EquilibriumVolume;
			return double.NaN;
		}
	}
}
